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ABSTRACT 

The metal budgets in some white dwarf (WD) atmospheres reveal that volatile-rich 
circumstellar bodies must both exist in extrasolar systems and survive the giant branch 
phases of stellar evolution. The resulting behaviour of these active asteroids or comets 
which orbit WDs is not well-understood, but may be be strongly influenced by sub¬ 
limation due to stellar radiation. Here we develop a model, generally applicable to 
any extrasolar system with a main sequence or WD star, that traces sublimation- 
induced orbital element changes in approximately km-sized extrasolar minor planets 
and comets traveling within hundreds of an. We derive evolution equations on orbital 
timescales and for arbitrarily steep power-law sublimation dependencies on distance, 
and place our model in a Solar system context. We also demonstrate the importance of 
coupling sublimation and general relativity, and the orbital consequences of outgassing 
in arbitrary directions. We prove that nongravitational accelerations alone cannot re¬ 
sult in orbit crossing with the WD disruption radius, but may shrink or expand the 
orbit by up to several au after a single pericentre passage, potentially affecting subse¬ 
quent interactions with remnant debris and planets. Our analysis suggests that extant 
planets must exist in polluted WD systems. 

Key words: minor planets, asteroids: general - stars: white dwarfs - methods: nu¬ 
merical - celestial mechanics - planet and satellites: dynamical evolution and stability 
- protoplanetary discs 


1 INTRODUCTION 

Water is known to exist in Solar system planets, moons, 
comets and asteroids. Even after almost 5 Gyr of Solar evo¬ 
lution, active dynamical processes inv olving water, like th e 
ventilation of plumes on Enceladus (iHansen et al.l l2006ll . 
persist. Comets still routinely shed their volatile content 
during close approa ches to the Sun. T he recent discovery of 
active asteroids fsee ljewitt et al.ll2015l . for a review) demon¬ 
strate that significant amounts of ice can exist in these bod¬ 
ies. But what will happen to the ice and water when the Sun 
turns off of the main sequence? 

The giant branch (GB) phases of evolution are violent, 
featuring stellar mass loss and intense radiation. In principle, 
this radiation will evaporate any surface water content of or¬ 
biting bodies out to several au, a l thoug h internal volatiles 
may be retained (lJura fc XrJl201Cll . 12012 1. In this paper, we 
refer to these orbiting bodies as two types: comets and minor 
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planets. In the context of exoplanetary systems, we define 
comets as bodies that formed beyond the snow/frost line, 
the line at which water and other volatiles froze out during 
system formation. The volatile content of comets may ap¬ 
proach 100 per cent, and represent a complex cocktail of e.g. 
water ice, carbon dioxide, carbon monoxide, and methane. 
Alternatively, minor planets formed closer to the star, and 
hence are less volatile-rich. 


Whatever events do occur during the GB phases, we 
have concrete observational evidence of the final outcome af¬ 
ter the stars have become white dwarfs (WDs): atmospheric 
meta l pollution w i th a high volatile content (iFarihi et al.l 
I2OI3I : iRaddi et al. I2OI5II. In the case of the WD named 
GD 61, Farihi et al.l ~( 201^ found that the progenitor of the 
debris in th e atmo sphere was 26 per cent water by mass; 
iRaddi et al.l (l2015l f later found that the atmosphere of the 
WD named SDSS J124231.07-1-522626.6 is enriched with el¬ 
ements that arose from a progenitor that was 38 per cent 
water by mass. 

Between one quarter and one half of all WDs harbour 
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metal-polluted atm ospheres (IZuckerman et aP l2003l . I 2 OI 0 I : 
iKoester et al.ll2014h . The dense natur e of WDs and the re- 
sulting fast sinking times (e.g. Fig. 1 of lWvatt et akllioidh in 
their atmospheres indicate that the currently-observed met¬ 
als do not arise from GB remnants nor internal dredge-up. 
Further, the spatial distribution and density of interstellar 
medium c louds prevent the metal pollution from originat¬ 
ing t here llAannestad et al.|j 19931: FViedrich_et al.l [2004l : ljural 
l200fil : [Kilic fc Redfieldll2007l : lFarilii et al.ll2nifll l. Instead, the 

metal pollution must arise from circumstellar material; bol¬ 
stering this hypothesis is the presence of orbiting dusty and 
gaseous discs. 

All 37 discs discovered around WDs so far contain dust 


lIZuckerman & BecklinI 19871: Becklin et al.l 20051: Kilic et al. 

20051: 

leach et al.l l2005l: iFarihi et al.l 20091: Berefors et al. 

20 14 

Rocchetto et al.l I 2 OI 5 II: 7 of them also contain 


strone easeous components (IGansicke et alJ 20061. 

2007|, 

20081: iGansic 

cell201ll: iFarihi et al.l l2012l: iMelis et al. 

2 OI 2 I: 

Wilson et al.l 

2 OI 4 I: Manser et al.ll2015l: IWilson et al.l 

2 OI 5 I). 


These discs exist in exceptionally compact conhgurations 
(~ 0.6 — 1.27?0, where Rq is the Solar radius) that are un¬ 
like anything which is observed in the Solar system, and 
are tho ught to arise from t he tidal d i sruption of minor 


planets ([Grah nnietalJ Ju^ 20^; IPebes et al.l I2OI2I : 


iBear fc Soketl 20131 : Veras et al. 2014ali and subsequen t cir¬ 


cularization of the resulting debris ( Veras et al.ll2015^ '). The 


disruption of comets is assumed to provide a smaller con¬ 
tribution to these discs, but is important otherwise for de- 
livering large quantities of hydrogen to WD atmospheres 
dVeras et aI.ll201^ : IStone et al.ll2015l ~). 

The minor planets or comets disrupt by entering the 
disruption (or Roche) radius - located at a distance of ap¬ 
proximately li ?0 - from a relatively great distance of at 
least several an, and hence contain eccentricities near unity. 
The semimajor axes of the minor planets or comets must be 
several an because otherwise they would have been engulfed 
by the expanding stellar envelope of the GB star. Surviv¬ 
ing grains, pebbles, comets, and minor and major planets 
are affected by a variety of other for ces which may per - 
turb these bodies to greater distances (I Veras et al.l 1201 ^ 1 . 
Planets are usually to o large to be affected by forces other 


than stellar mass l oss dOmarov 


IVeras et al.l I2OIII : lAdams et al 


19621: Hadiidemetriou 196^: 

l2013nVeras et al.l l2013al l. 


Minor planets and comets, how ever, will also be per¬ 
turbed by th e Yarkovsky effect dVeras et al.l l2015al l. the 
YORP effect dVeras et aL l201 4bll. Povnting-Robertson drag 
and ra diation pressure ( Bonsor fc WvattI 2010|;_ Veraset_^ 
2015al~l. and possib ly stellar wind drag ( Dong et al. 2010l : 


ions iKratter & PeretsI I 2 OI 2 I: IVeras & Tout! 2012lf 

Dr sur- 

viving planets in the system ( 

Debes & Sieurdsson 

2002: 

Bonsor et al.l I 2 OIII: iDebes et al. 

2012: 

Veras et al. 

2013b: 

Vovatzis et al.l 20131: Frewen & HansenI 

2014: Mustill et al. 

2 OI 4 I: Veras & Gansicke| 2015|'). 


Such complications ensure that minor planets and 
comets will exist on a wide variety of orbits around WDs, 
and not necessarily ones which enter the WD’s disruption 
sphere. So what happens to the volatile-rich minor planets 
and comets which miss the disruption sphere? In this work, 
we isolate the orbital effects of the sublimation and out- 
gassing of volatiles due to WD radiation, and focus on the 
minor planets with highly eccentric orbits that fail to reach 


the WD disruption radius. Our results are also applicable 
to a broad variety of comets, as long as their volatile con¬ 
tent is not high enough to lead to break up or non-periodic 
behaviour. We create orbital models which are widely ap¬ 
plicable to nearly any extrasolar system assuming that no 
major planets interfere with the minor planet’s orbit. 

In Section 2, we outline our model assumptions for 
volatile sublimation. Sections 3 contains the resulting un¬ 
averaged orbital element evolution equations, along with a 
few applications, demonstrating both the interplay of gen¬ 
eral relativity and sublimation and the nearly-fixed nature 
of the orbital pericentre. We then reconsider how this peri- 
centre might change due to outgassing in Section 4, where 
we also relate the extrasolar and Solar system sublimation 
models. We discuss the implications in Section 5 and con¬ 
clude in Section 6. The Appendix contains an extension of 
our sublimation model that removes the key assumption of 
a specific power-law exponent. 


2 SUBLIMATION MODEL BASICS AND 
MOTIVATION 


Genturies of observations of Solar system comets reveal or¬ 
bital changes at each apparition. These changes, which can 
represent shifts in semimajor axis and perihelion of tenths of 
an au on short time scales (e.g. see Section 5 and Figure 3 of 
ISzutowi^l200Cll . and lMaauetlboi^ may have dramatic long¬ 


term consequ ences. Sungrazing comets often are destroyed 
outright (e.g. Brown et al.ll201ll ~). such as G/2012 SI ISON 
dKnight fc Battamsll201 J T Orbital alterations are primarily 
caused by two factors: (1) gravitational interactions due to 
the architecture of the Solar system, which contains many 
major and minor planets and moons, and (2) nongravita- 
tional forces. 

The presence of volatiles ensures that the nongravita- 
tional forces from sublimation and outgassing will dominate 
those from radiation pressure, Poynting-Robertson drag, the 
Yarkovsky effect and YORP. By sublimation we refer to the 
mass exodus of surface particles which are distributed ap¬ 
proximately evenly over large parts of the minor planet or 
cometary surface. By outgassing we refer to localized vio¬ 
lent eruptions of material in random directions. This sec¬ 
tion and Section 3 covers sublimation, and Section 4 covers 
outgassing. 

Minor planets which are strongly affected by sublima- 
tion or outgassin g are known as “active asteroids” (e.g. 
Ijewitt et al.ll2015II . Glassically, strong outgassing as seen in 
some comets takes the form of a visibly dramatic tail. Mod¬ 
ern observations reveal that the situation is much more com¬ 
plex, s ometimes including up to six tails for an individual 
object (|jewitt et al.ll201^ . 


Consequently, small-body sublimation and outgassing 
has predominately been investigated observationally in the 
context of the Solar system and from the perspective of an 
observer on Earth. Theoretical developments have largely 
followed suit. Here, however, we consider sublimation from 
a more general perspective, for application to exoplanetary 
systems. For WD systems, we have not yet directly detected 
comets, major planet^, moons nor any intact minor planets. 


^ The object WD 0806-661 b is by many definitions a bona- 
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Sublimation of planetesimals in white dwarf systems 3 


and hence cannot utilise many of the tools developed for or 
assumptions about the Solar system. 


2.1 Physical properties of the minor planet 

Denote the mass of the minor planet as M, which contains 
both a rocky non-volatile component of mass Me and a 
volatile component of mass Mv such that 


M = Mc + Mv. (1) 

The non-volatile content may be thought of as the minor 
planet’s core, although we need not make a priori assump¬ 
tions about just where the volatiles are located. We consider 
only minor planets that will remain largely intact, such that 
they do not fragment or sublimate away a large fraction of 
their mass. Consequently, the assumption that we do make is 
that Me Mv, so that M « Me remains constant over the 
number of pericentre passages needed for all of Mv to subli¬ 
mate away. This assumption is robust for minor planets and 
many types of comets. Alternately, for some comets, Mv/M 
may approach unity. In this case, the dynamical evolution 
can change in more drastic ways than we report here. 

The rotation state of the minor planet will crucially 
determine where outgassing will occu r, with consequences 
for the minor planet’s orbital evolution. ICheslev fc Yeomans! 
ll2005ll have shown that outgassing via jets can be modeled 
in an average sense, if the rotation of the object is well de¬ 
fined. Their model suggests that outgassing is predominately 
happening on the starlit hemisphere on slow rotators with 
low thermal inertia. Yet, outgassing it self will likely chang e 
the rotation state of the minor planet llMottola et al.ll20l4) , 
thereby randomizing the perturbations and emphasizing the 
stochasticity of the process of outgassing. Stochastic varia¬ 
tions intrinsically contain time-dependent variations in di¬ 
rection and magnitude, as would occur in a tumbling minor 
planet. The a ctual accelera tions will be probabilistically dis¬ 
tributed (e.g. |Pierredf2014l ~). 

Treating potentially complex rotation is beyond the 
scope of the present work and beyond the scope of our 
knowledge of minor planets in extrasolar systems. Instead, 
we consider outgassing from a general perspective in Section 
4, rather than ass uming outgassing oc curs only on the starlit 
hemisphere, as in iFestou et al.l lll993ll . 


2.2 Orbital change of the minor planet 

Many forces contribute to the orbital motion of the minor 
planet, including those due to radiation, sublimation and 
outgassing. Further, minor planets with orbital pericentres 
close to the WD might be significantly affected by general 
relativity (GR). If the minor planet has an orbital period 
similar to those of exo-Oort cloud comets, then near apoc- 
entre both Galactic tides and stellar flybys can cause signif¬ 
icant perturbations. 

Strictly, when a body (like a minor planet) loses mass 


fide planet orbiting a WD at a distance of about 2500 au 
llLuhman et al.ll2011 i. Further, an increasingly strong case is be¬ 
ing made fo r a pair of purport ed planets which orbit a WD - M 
star binary llMarsh et al.ll2oI3l . 


anisotropically, its orbit is perturbed by both linear momen¬ 
tum r ecoil and mass loss . Although fully described analyti¬ 
cally llVeras et al.ll2013D h the equations of motion due solely 
to mass loss can be ignored here because M <C Mwd, where 
Mwd is the mass of the WD. The overall change in mass of 
the system due to sublimation plus outgassing is negligible. 
The primary perturbations on orbits due to sublimation and 
out gassing is li near momentum recoil. 

Whipplj lll95df used conservation of momentum to 


show that the nongravitational acceleration on a minor 
planet (Jng) is related to its ejecta through 


MJng = —DrrisVs, (2) 

where ms is the mass of one molecule of the molecular species 
that is sublimating, Vg is the mean ejection velocity of that 
species, and D is the number of molecules of that species 
sublimating per time. For our numerical computations, we 
will assume the sublimating species is water ice (H 2 O), and 
hence will adopt rris = 2.99 x kg. We assume that 

most of the sublimation occurs on the starlit hemisphere 
because the minor planet’s tiny ther mal inertia will ensur e 
a near-zero thermal wave lag angle (iGroussin et al.ll2013l l. 
Gonsequently, the complex physics which relates thermal 
balance, rotation, sublimation and outgassing is simplified 
here through the values of D and Vg. 

The equation of motion for the minor planet is 


G (Mwd + Me) f 


r^sub “t“ r extra, 


( 3 ) 


where 


r'sub '— 0 (A/v) 


( 4 ) 


r is the distance between the minor planet and the WD, 
^sub is the acceleration due to sublimation, Fextra is the 
acceleration due to other forces, and 0 is the Heaviside 
step function such that 0(Mv) = 1 when Mv > 0, and 
zero otherwise. Other forces which must be included in 
Fextra depend on the system being studied. These additive 
terms from other forces have been used or applied in previ¬ 
ous investigations: mass loss dOmaro ^ Il962l : riladj jdemetrjoul 
Veras et aI]l2013a^ . the Yarkovsky e f fect dVeras et al ' 


2015a l. t he YORP effect dVer as et al.| |2014b |), Galac¬ 
tic t ides (iHeisleLjMHCTnahlJ 198^ RFdomn ^^jVoeschld 


1 19871 : iFouchard et al.l I 2 OO 6 I : IVeras fc Evans! l2013bf l and 
stellar flybys in the Galactic field (iFreeeau et al.l l2004l : 
IZakamska fc Tremainel [20041 : IVeras &: Moeckel 2012l l. 

One effect which might strongly influence the orbit of 
any minor planet or comet with a pericentre close to the star 
is GR. Solar system studies have indicated the importance 
of GR in calculations of cometary orb ital solutions (e.g. 
lYeomans et al.iri996l : iMaauet et al.ll2012ll . Consequently, we 
do incorporate GR into our formalism. We do not consider 
the other effects mentioned above any further, and hence set 
Fextra = Fgr. This assumptiou is reasonable because large 
minor planets with low thermal inertia are not strongly af¬ 
fected by Yarkovsky forces. Similarly, sublimation and out¬ 
gassing dominate the evolution of minor planet spin as long 
as volatiles are present. 

We use the complete expression for the IPN (leading 
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4 Veras, Eggl & Gdnsicke 


post-Newtonian) term of GR i n the limit o f a zero-mass 
minor planet (equation 3.190 of lBeutl^l2005h 


rcR 


GMwd 

q 2 j .3 


^ 4GMwd 



r + 4 (u • r) n 


(5) 


where v is the velocity of the minor planet with respect to 
the star and c is the speed of light. 

Now consider equations equations in greater 

depth. In order to develop physical intuition for how the 
2-body orbit changes due to linear momentum recoil, we in¬ 
sert equation into equation and reexpress equation 
in terms of orbital elements. Doing so requires us to ex¬ 
press Jng in terms of positions an d velocities. In general, 
D oc r~^ (ISosa fc Fernand^ l2009l ~). as the gas production 
rate peaks at the pe ricentre. More complex expressions (e.g. 
iKorsun et al.l l2014h are Solar system-specific, and are de¬ 
pendent on_observables. Also, Vg oc (equation 28 of 

IWhippla I 1950 I I. We note that this velocity dependence is 
differ ent from the dependence predicted by equation 

8 of llKorsun et al. Ilioia) because the latter is the terminal 
velocity of the gas, rather than the velocity at the area of 
the ejection. Assuming that sublimation occurs only on the 
starlit side, we find Vg/vg = rjr. 

These dependencies yield: 


main conclusion without forcing us to resort to specific pa¬ 
rameter choices. 

Also, for notational ease, we will henceforth drop the 
0(Mv) term in all subsequent equations. Nevertheless, we 
must determine over what timescale Mv is depleted. The 
volatile budget decreases according to 


Mv = M^{0) - m^Do dt 


Mv(0) - 


rUaDorQ 

— e 


(/-/o) 


(7) 


Mv(0) - 8.6 X 10 '^kg X N 0-999 7 ^ 

/ a \ - 1/2 / ro Y f M \ 
llOau/ ll0“^au/ \0.6 Mq/ 

Do 


- 1/2 


2.99 X 10~^®kgy \10^®mol/s 


(8) 


where N is the approximate number of orbits, n is the mean 
motion, a is the semimajor axis, e is the eccentricity and / 
is the true anomaly. The number of orbits over which the 
entire volatile budget would be depleted is 


rsub = -0 (Mv) 


DoVgoiris 

Me 


r 

\ r J r 


( 6 ) 


where D = {ro/r)^Do and Vg = Ugo(ro/r) ^ (r/r). Although 
these power-law dependencies are physically-motivated, we 
realize that other dependenc ies may be used bas ed on So¬ 
lar system observations (e.g. ISanzovo et alll200l| j. In order 
to take into account other possibilities, we generalise equa¬ 
tion (O in the Appendix. There we show (in Table [XT]) our 
qualitative results hold for steeper power-law dependencies. 

We cannot assume that the composition and outgassing 
properties of extrasolar minor planets are similar to Solar 
system minor planets because their formation channels may 
be different and because WDs and the Sun have different 
power spectra. Nevertheless, we can look towards the So¬ 
lar system for rough guidance. In the Solar system, rep¬ 
resentative Do values for sublimation are a few orders of 
magnitude less than than the observed outgassing-based val¬ 
ues of Do ~ 10^® — 10®^ molecul es per second (Tab le 3 of 
lA’Hearn et al.|[i995l and Fig. 3 of ISosa fc Fernandejl20ll]l . 
The ejection speeds |ug| from the surface are unknown, but 
the terminal escape speed of t he gas is thought t o be on 
the order of cm/s - tens of m/s (iKorsun et akll^ld l. If ro is 
taken to be at the pericentre of the orbit, then equation m 
effectively models the diminution the escaping molecules as 
the minor planet flies away from the WD. This sublimation 
model is adequate for our purposes. 

Real outgassing naturally does not occur in such 
a well-behaved radially symmetric manner; nonzero tan- 
gential components may be signiheant (e.g. T able 1 of 
iMarsden et al.lfi973l . Table 2 of Krolikowskall2003 . and Table 
4 of Szutowicz &: RickmarJ[2006l b Therefore, the initial con¬ 
ditions for exosystems with minor planets and comets are 
largely unconstrained by observationally-motivated consid¬ 
erations. Consequently, in Section 4 we consider outgassing 
from a general analytical perspective that strengthens our 


Mv (0) no ao ^/l - eg 
27rmsDorg 


(9) 


Equation ® is useful for determining simulation durations, 
and how much of the body’s mass has been lost when the 
body’s orbit begins to undergo major changes. However, that 
estimate becomes progressively worse as the semimajor axis 
undergoes drastic variations. 

Having obtained an observationally-motivated func¬ 
tional form for sublimation (equation [6]), we can now derive 
the equations of motion in orbital elements to understand 
how sublimation changes the orbit. 


3 UNAVERAGED EQUATIONS OF MOTION 
FOR SUBLIMATION 

As a first step, we consider the equations of motion with as 
few assumptions as possible. In the subsequent analysis, the 
following auxiliary set of variables will be useful: 


Cl = ecosuj + cos{f + uj), (10) 

C 2 = esincu-|-sin (/-I-oj), (11) 


Gs E 

-= (3 

-b 

4e 

COS 

/ + 

cos 

2/) sincu 

-b 2 (e -b cos /) 

cos UJ sin /, 










(12) 

Ge = 

= (3 

-b 

4e 

cos 

/ + 

cos 

2/) cos a; 

— 2 (e -b cos /) 

sin UJ sin /, 










(13) 

Cr = 

= (3 

-b 

2e 

cos 

/- 

cos 

2/) cos a; 

-b sin cu sin 2/, 

(14) 

Cs = 

= (3 

- 

cos 2/) sin 

iJ - 

- 2 (e -b cos /) cos UJ sin /, 

(15) 

Co = 

= (3 

-b 

2e 

cos 

/- 

cos 

2/) sincu 

— cos to sin 2/. 

(16) 


This nomenclatu re was chosen to main ta in consistency with 
previous studies (IVeras fc Evansl2013al f9: IVeras et al.ll2013^ . 


l2014dl : IVerasI l2014a ) and highlight the naturally-occurring 


© XXXX RAS, MNRAS 000.ITIII41 




















































Sublimation of planetesimals in white dwarf systems 5 


quantities in many formulations of the perturbed two-body 
problem, including those with sublimation. 


3.1 General orbital evolution equations 

In order to obtain the unaveraged equations of motion 
in orbital elements, we employ the proce dure suggested 
in IVer as et al.l ([2oT3) and embellished in IVeras fc Evan^ 

e ), b ased o n pert urbative analyses from lEfroimskv 
and lGur^ (l2007l l. We find 


a\ _ 2 DoVgomBrl^^ [I + e cos 

I/sub Mcna9/4 (1 - 

X Cl (cos i cos fl — cos i sin Q. -\- sin i) 


—C 2 (cos n -I- sin fl) 


(17) 


DoVgomsrl^'^ (1 -I- e cos 

2Munai3/4 (1 - 


X Ce (cos i cos 12 — cos i sin 12 -|- sin i) 

—C 5 (cos 12 + sin 12) 

di\ _ DoVgomarg^'^ {1 + ecos cos {f + Lu) 
‘t-^/suh Mcnai3/4 (1 _ e2)’^/4 

X sin i (sin 12 — cos 12 -|- cot i ), 


(18) 


(19) 


/ dl2 \ DoVgorUarQ^'^ (1 -I- e cos Z)®^"^ sin (/ -|- uj) 

/ sub Mcna^^/^ (1 — 


centricities near unity, sublimation-induced changes to the 
semimajor axis, eccentricity and pericentre distance are rela¬ 
tively strong, minor and negligible, respectively. Altering the 
pericentre distance of an orbit of a highly eccentric minor 
planet around a WD through sublimation alone is difficult. 
Other agents in the system are needed to perturb the minor 
planet into an orbit where it can tidally disrupt. 

We can facilitate the coupling of GR with sublimation 
by considering the corresponding equations of motion in or¬ 
bital elements. By ass uming M <C Mwd to an excellent 
approximation we find (IVerasll2014bh 


f da\ ^ 2G^MwDesin/(1-|-ecos/)^ 
V'^^/GR c^na^ (1 — 

X [7 -|- 3e^ -I- lOe cos /] , 

/de\ ^ G^Mwd sin / (1 -I-ecos /)^ 
V'^^/GR c^na"^ (1 — 

X [3 -I- 7e^ -I- lOe cos /] , 



(24) 


(25) 

(26) 


/ duj\ ^ G^(1-t-ecos/) 

\ di / gr c^na'^e (1 — e^)®^^ 

X [(e^ — 3 ) cos f + Se — 5e cos (2/)] , (27) 

/dqX _ G^Mwd sin/(1-I-ecos/)^ 

\d2/GR c^na^ (1 — e)®^^ (1-I-e)^^^ 

X [-3 -b 8e -b - lOe cos /] . (28) 


X (sin 12 — cos 12 -b cot i ), 


( 20 ) 


diu\ _ DoUgomsrQ^"* (1-b ecos Z)®^"^ 

d-i J sub 2 Mcnai 3 / 4 g 

X Gs (cos i cos 12 — cos i sin 12) -b Gg sin i 
-bGr (cos 12 -b sin 12) -b 2e sin (f + uj) cos i cot i 


( 21 ) 


f dq\ _ DoUgomarQ^"* (1-becosZ)®'^'^ . f f\ 

Vdi Aub ^ Muna9/4 (1 - e)®/^ (1 + e)®®/^ VV 


3.2 Planar orbital evolution equations 

Equations ([T7ll-([22ll demonstrate how sublimation affects 
the entire set of orbital elements, including the inclination 
and longitude of ascending node. However, because WDs are 
extremely spherical compared to planets, the orbital archi¬ 
tecture is rotationally symmetric. Hence, future WD studies 
might need to consider just the planar case. In order to ob¬ 
tain the equations for the planar case, we take the limit 
of zero inclination, and define the longitude of pericentre 
zj — Q + ui. 


— 2 cos ( ^ ) (2 -b e — cos f) 


cos UJ (cos 12 -b sin 12) 


-b sin UJ (sin i + cos i (cos 12 — sin 12)) 
sin UJ (cos 12 -b sin 12) 
— cos UJ (sin i — cos i (sin 12 — cos 12)) 


+4 sin® f I 


For the planar equations of motion, the G variables 
which contain a subscript of “P” (below) are equivalent to 
those variables dehned in equations (HSli-dleJ except that uj 
is replaced with vj. With these definitions, we can compactly 
express the coplanar equations of motion as 


"I 'I 

/ da\ 

(22) 

[d^J 


p 

sub 


2 DoVgomBrl^^ (1 -b e cos 

Mc'na^/'^ (1 — 


(Gip — Ggp), 


df n(l-becosZ)^ duj .dil 
dt (1 _ e^)®'^^ dt dt 


(23) 


where i is the inclination, 12 is the longitude of ascending 
node, UJ is the argument of pericentre and q is the pericentre 
distance. Note that the form of equation (1231) . without the 
subscript “sub”, is not restricted to sublimative forces. 

Our main result can be hinted at directly by inspection 
from equations (II3-(I23; for minor planets with orbital ec- 


(29) 



DoUgornsrp^"* (1 -b e cos 
2Mcnai3/4 (1 - 


(Gep — Gsp), 


(30) 



DoVgomsry'^ (1 -b e cos Z)*^"^ 

2 Mcnea^^/'^ (1 — 


(Ctp -b Ggp) , 


© XXXX RAS, MNRAS OOO.fTHTil 
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(31) 


_ DoVgomsr®^'^ (1 +ecos/)®^'‘ 

dt) 2Mcna9/4 (1 - ef'* (1 + 


X 


cos (2/ + w) + cos vj (3 


4 cos / + 2e sin /) 


— (3 — 2 (2 + e) sin /) sin-1x7+4 sin (/ + w)—sin (2/ + zo) . 

(32) 


Equations ([^-(|32l) show that even in the planar case, 
the orbital equations are still nontrivial functions of the true 
anomaly of the asteroid. Nevertheless, we can glean some 
physical insight from the formulae. The eccentricity terms in 
the denominator illustrate that {da/dt)sub will evolve more 
quickly than {de/dt)s^h, {dvj/dt)suh, and {dq/dt)sub- 

For GR, because equations (|24l)-([25l) and (l^ - (l28l) are 
all independent of i, and w, we have 


( 

( 

( 

( 


da 5 

P 

dt ) 

GR 

de\ 

P 

dt J 

GR 


P 

dt J 

GR 

dzu 

Y 

dt 

/ GR 



(33) 

(34) 

(35) 

(36) 


3.3 Numerical integrations 

In this section we provide some specific examples of minor 
planet orbital evolution due to sublimation and GR. 


3.3.1 Setup 

First, we must establish some fiducial parameters. A crucial 
parameter is Dq , which varies by several orders of magnitude 
in the Solar system. This variation should be even higher in 
WD systems because WD luminosities and spectral prop¬ 
erties vary by orders of magnitude depending on the WD 
cooling age; WDs may be either more luminous or (much) 
less luminous than the Sun. Although the latter case holds in 
the majority of cases, Dq cannot necessarily also be assumed 
to be lower because minor planets can reach much smaller 
values of q in WD systems. In fact, without being engulfed, 
a minor planet may attain q = 6 x 10“® au. However, this 
value is well within the WD disruption radius, where the 
minor planet would likely be shorn apart. Instead, q = 0.02 
au would nearly guarantee that the minor planet survives; 
actually survival is likely to be ensured over hundreds of 
orbits for val ues of q which are one or two orders of mag¬ 
nitude lower JVeras et al.l 120143) . Gomplicating these con¬ 
siderations for determining a fiducial Dq is that sublimation 
is a function of temperature and pressure, and significantly 
localized physical properties such as shadowing, thermal in¬ 
homogeneities, and spin. 

Gonsequently, we defer linking Dq with detailed in¬ 
ternal minor planet models and WD cooling age to other 
studies. Instead, here we identify the value of Dq in which 


sublimative-induced orbital perturbations becomes signifi¬ 
cant for a minor planet which does not initially quite reach 
the WD disruption radius (for example, qo — 0.02 au). Like 
Dq, the value of Ugo is intrinsically linked to local proper¬ 
ties of the minor planet and is unknown outside of the Solar 
system. Fortunately, we can partially circumvent our igno¬ 
rance by recognizing that all of the sublimation-based orbital 
evolution equations are directly proportional to the product 
DoVgQ. This property allows us to reduce the number of de¬ 
grees of freedom in any phase space exploration. For our 
explorations, we choose the four products corresponding to 
the range DoVgo = mol/s) x (Ikm/s). These speeds 

are easily achieved; in the Solar system, assuming a Maxwell- 
Boltzmann distribution, one obtains (ugo) = 675 ± 275 m/s 
(with e xtreme speeds of sev eral km/s) for temperatures of 
360 K llGroussin et ahllioistl . 

Our other adopted val ues include Mwd = Q.QMp, (rep¬ 
resenting a typical WD; see|Ljebert_et_al7 [200d; iGamenzindl 
I 2 OO 7 I : iFalcon et al.l l2010l : iTremblav et al. 2013li , and Me = 
10^® kg (roughly corresponding to an active asteroid with a 
radius of 1 km and density of 2 g/cm®). We place our mi¬ 
nor planets on ag = 10 au orbits, because most km-sized 
minor planets within about 7 au wi ll have been spun up to 
fission from giant branch radiation (|Veras et al.ll201^ L We 
do adopt q — 0.02 au, meaning that eo = 0.998, a typical ec¬ 
centricity for a minor planet just missing the WD disruption 
radius. We sample minor planets on both slightly inclined 
(~ 1°; Figure [T]) and highly inclined (~ 45°; Figure [2| or¬ 
bits with respect to the same fixed but arbitrary reference 
plane. Our simulations begin at the pericentre and we prop¬ 
agate the systems for approximately 25 orbits. We integrate 
equations HZll-dlT]) and (I23l)-([27l) all together. 


3.3.2 Results 

We present our results in Figs.[T][21 The sawtooth-like curves 
in both figures indicate that sublimation causes the greatest 
orbital changes at the pericentre, so that each step corre¬ 
sponds to one orbit. 

Figure [T] showcases a side-by-side comparison when GR 
is not included (left panels) and is included (right panels) in 
the simulations. The difference is striking, particularly for 
the semimajor axis and pericentre evolution; for the blue, 
pink and green curves, at least, GR dominates over subli¬ 
mation. Alternatively, the black curve for the evolution of 
07 indicates a regime where sublimation may dominate GR 
(after a few hundred yr). The figure makes clear that for 
reasonable parameters that one would expect in a WD sys¬ 
tem with planetary remnants, GR can play a crucial role 
and should be included in sublimation models. 

Another comparison to consider in that figure is be¬ 
tween the top and middle plots, whose j/-axes are equivalent 
and in au. The change in q — qo is typically orders of magni¬ 
tude less than the change in a — oq, both with and without 
GR. GR dominates the evolution of q in most cases, flat¬ 
tening out the curves. In no case does q vary by more than 
10”'^ au. 

In Figure[^ we include GR and show the time evolution 
for e, i, and D and as well as for a, q and u, but now for 
minor planets with different initial values of i, and uj. The 
highly inclined retrograde orbits sampled here yield different 
evolutions from Figure 1, but reinforce the main result: the 
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Without GR With GR 





Figure 1. The orbital changes over about 25 orbits experienced by a 10^® kg minor planet orbiting a 0.6Mq star (like a WD) and 
sublimating water ice at a given rate Dq and with a given speed Ugo- The product of Dq and Vgo scales linearly with the change in 
each orbital element. The left panel does not include GR, and the right panel does. We set the initial minor planet orbit at ao = 10 au, 
qo = 0.02 au, eo = 0.998, iq = 1.0°, Qq = 180°, uiq = 1.0° and /o = 0.0°. The plot demonstrates the importance of including GR when 
the pericentre is close to the star, and hints that q cannot drift into the star’s disruption radius by the combined forces from sublimation 
and GR alone. 


pericentre change does not exceed 2 x 10“"^ au for any curve, 
whereas the semimajor axis may vary by as much as 0.05 au. 
Note that the two uppermost plots are drawn on the same 
scale. 


For both sets of simulations, we also computed the to¬ 
tal amount of mass actually lost by the minor planet from 
equation © throughout the evolution. This value is about 
9.7 X 10® kg, four orders of magnitude less than than the 
core mass of the minor planet. 


4 TRANSVERSE AND NONPLANAR 
OUTGASSING 

So far our model for sublimation assumes that the nongrav- 
itational acceleration acts opposite to the radial direction. 
But what if the minor planet experiences a significant out- 
gassing episode in another direction? In order to answer this 
question, we first consider a commonly-used formulation for 
outgassing in the Solar system. 


4.1 Sublimation in the Solar system 

This prescription effectively models a perturbed two- 
body problem, where the perturbative force is the prod¬ 
uct of a semiempirical functional form for water ice 
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8 Veras, Eggl & Gdnsicke 


With GR 


f/yr 



t/yr 



t/yr 


Figure 2. Like Figure [T] except with io = 45°, Ho = 200° and wq = 135°. Both sublimation and GR are included in the simulation. 
Illustrated here are the time evolution of all of the orbital elements except /. Despite the minor planets having highly inclined orbits, q 
remains robust to major orbital changes. 


(|DeIsemme fc Milleil Il97ll l and a linear combination of 
the r adial, transverse an d normal comp onents of the mo¬ 
tion llMarsden et al.lll97^'l. Equation 1 ofISosa fc Fernanded 
l| 2011 fl and equation 1 of lSzutowicd (| 200 nl l demonstrate that 
the Marsden formulation of outgassing for Solar system 
comets may be expressed as 

G {Mq + Me) 


+N(r) 


N(r) = a 




+^2 ■ 


r X II 


Radial 
= R 


ro 


1 + 


■ +A3 


r X V 


Transversal 

= f 



(37) 


(38) 


such that the constants 77 , ^ and ^ are empirically derived, 
with typical fiducial values of 2.15, 5.093 and 4.6142, respec¬ 
tively. Further, the values of Ai , A 2 and A 3 are not constant 
with time; for real comets, they can vary by an order of mag- 
nitude or more (see Table 2 of Froesch le fc Rickman|[l9^ 
and Table 2 of ICheslev fc Yeomans! 2005h . Even this for¬ 
mulation, however, cannot explain many observational phe¬ 
nomena. For example, sometimes outga ssing is not symmet - 
ric about the pericentre. Equation 5 of iKrolikowskal (l2004h 
takes into account such cases with related functional forms. 


Note that we recover the formulation used in Sec¬ 
tion 3 by setting A 2 = A 3 = 0, r/ = 9/4, (/ = 0 and 
a = 0 (Mv) DoVgonis/Mc. We find the equation for the semi¬ 
major axis variation from the Marsden formulation to be 



^ 2 [^iesin/-b A 2 (l-becos /)] 

n(l-e2)^/^ 


(39) 


out 
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Sublimation of planetesimals in white dwarf systems 9 


This form is equival ent to equation (24) o f Bunr^ lll976l i (and 
equation (2.145) of iMurrav fc Dermottlll999ll . but without 
the factor of N(r), even though we made no assumption 
about the magnitude of the nongravitational force (the force 
in those other studies is assumed to be small). The fac¬ 
tor of N(r) does not alter the classical functional form ob¬ 
tained due to the radial, transverse and normal components. 
The evolution of the other orbital parameters are as follows 
and c an also be re-expressed in the forms present in iBurn^ 
lll976tl . 


f) =N(r)x 




2na (1 -I- ecosf) 


2Ai sin / (1 -I- ecos /) -I- A 2 (4 cos / -I- e (3 -I- cos (2/))) 


rfi/out an (1 +ecosf) 


^ CSC»in(/ + a.)^ 

dt J an (1 +ecosf) 

dij\ 

dt J . na 

A out 

cos/ sin/ (2-I-ecos/) cotisin(/ + w) 

- ^2 - 7 :r“- — - r ^3 


(40) 

(41) 

(42) 


e (1 -I- ecos /) 


1 -I- e cos / 


(43) 


From equations (I39ll - (l4n|l . we derive the change in peri- 
centre g and apocentre Q as 

=N(r)x - - 

dt J aut 2nVl - e^ (1 -|- ecos /) 

—2Ai sin / (1 -I- e cos f)+A 2 (4 -I- e — 4 cos / — e cos (2/)) 

(44) 


dQ 

dt 


(1 e)®^^ 

= N(r) X - \ - 

out 2nVl-e(l-f ecos/) 


2Ai sin / (1 -|- ecos f)+A 2 (4 — e -I- 4cos / -|- ecos (2/)) 

(45) 


We see that da/dt,de/dt,dq/dt and dQ/dt are unaf¬ 
fected by the choice of A 3 . If A 3 = 0 then the inclination 
and the longitude of the ascending node stay constant. If 
the motion is fully co-planar, i.e. A 3 = 0 and i = 0, then the 
ascending node is no longer well-defined and the evolution 
of the longitude of pericentre becomes 


d'uj . \ 

— ) =N(r)x 




aen (1 -|- ecos /) 


— Ai cos / (1 -I- e cos /) + A 2 sin / (2 -I- e cos /) 


(46) 


The formulae in this section allow one to model the evolu¬ 
tion of a particular exosystem minor body for an assumed 
set of {Ai, A 2 , A 3 }. Currently, no such sets have been obser- 
vationally measured outside of the Solar system. 


4.2 Invariance of the pericentre 

Regardless, we can reveal important properties of the mo¬ 
tion by analyzing the equations alone. We know that orbital 
changes predominately occur during pericentre passages; 
consequently, we can expand the equations about / = 0 °. 
Further, in WD systems, any minor planets or comets which 
closely approach the WD must have extremely eccentric or¬ 
bits (e 0.998). Therefore, we can also Taylor expand the 
equations about e = 1 —e, where e <C 1. The following result 
shows the explicit dependencies on (1 — e), Ai, A 2 and A 3 . 
We facilitate the expressions by first defining 




ro 


1 -f 


a (1 — e) 


ro 


1 -c 


(47) 


We obtain 


f about 0° 
e about 0 


A7{-^(l-e)-5 A 2 +O [(l-e)^]} 

+^{-^ (1 - e)-5 Ai + O [(1 - e)^] I / + O [f] . 

(48) 


Hence, the semimajor axis is most strongly influenced 
by a nonzero A 2 term, at the (1 — e)~^^^ order. The next 
strongest influence comes from a nonzero Ai term at the 
/ (1 — e)~^^^ order. The eccentricity dependencies are simi¬ 
lar, except weaker by a factor of (1 — e): 


f de\ 



f about 0° 
e about 0 


AT 


^(l_e)5 A 2 -bO 

an 


(1 



+A^ (1 - e)i Ai + O [(1 - e)i] | / + O [f] . 

(49) 


The inclination and longitude of ascending node terms 
are independent of Ai and A 2 . Instead, these orbital pa¬ 
rameters vary due to a nonzero A 3 , at the (1 — e)^^^ leading 
order as follows 



out 


/ about 0° 
e about 0 


AT 


COSU} 

y/2an 


(1 - e )2 A 3 -bO 



+Afi^^^il-e)i A3 + 0[{l-e)^]^f + 0[f], 

(50) 



/ about 0° 
e about 0 
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10 Veras, Eggl & Gdnsicke 


Af 


V 2 an sin i 

COS UJ 

\/ 2 an si 


(l-e)^ ^3 + 0 [(l-e)5] 

.(l_e)i43+0[(l-e)i]|/+0 [f]. 

(51) 


Unlike the evolution of a, e, i, and fl, the evolution of 
UJ is instead dominated by the combination of an Ai term 
and a A 3 term, as follows 


duj \ 


Af 


f about 0° 
e about 0 

2 _ , 1 ^ cotisincj _ , i ^ 

(1 - e) 2 Ai H- r= -(1 - e) 2 ^3 


\/ 2 an 
+0 [(l-e)i] 


V2a 


./-] 3 .1 . coticosoj _ ,1 ^ 

+Af <1 -^— (1 - e) 2 ^2 H- -j= -(1 - e) 2 A3 


\pian \p2an 

+0 [(l-e)i] |/ + 0 [/ 2 ]. 


(52) 


The coplanar version is 


dzu 

dt 


f about 0° 
e about 0 


+0 [f]. 


(53) 


Finally, now we compare the above results with the evo¬ 
lution of the pericentre 


= 0 


f about 0° 
e about 0 




(54) 


Equation (I54II demonstrates that (i) the pericentre under¬ 
goes no evolution at all in the leading order in true anomaly 
(/°) term, (ii) in the next order term (/^), the eccentricity 
dependence is weaker [(1 — e)^^^] than that for any other 
orbital parameter, and (iii) this weaker leading-order term 
is independent of A 2 and A 3 , so that the best way to change 
the pericentre is through radial perturbations. The first 
two points illustrate how robust the orbital pericentre is to 
changes compared to any other orbital element. 


5 DISCUSSION 

Our conclusion is strong and independent of the detailed 
changes which a minor body or comet may experience after 
each pericentre passage if the lost volatiles do not consti¬ 
tute a large fraction of its total mass. In this case neither 


outgassing nor sublimation can perturb a WD exosystem 
minor planet or comet into the WD disruption sphere if the 
original orbit remains outside of that sphere. 

Although potential time dependencies of many of the 
parameters that we held constant, including Dq and Ugo, 
could certainly change a minor planet orbit, the actual peri¬ 
centre cannot drift inward enough to transform an other¬ 
wise “safe” orbit into a “disruptable” one. Consequently, 
the presence of planets is necessary to perturb smaller bod¬ 
ies that reside within a few hundred an of the WD inside 
of the WD disruption sphere. In main sequence exosys¬ 
tems, the invariance of the pericentre distance also holds 
as long as the extrasolar minor planet or comet or is on 
a highly eccentric orbit. These systems do not necessarily 
need to host currently observable planets, as previous radial 
incursions from planet-planet scattering events can gener¬ 
ate h igh eccentricities of sma l ler bodies that survive ejec¬ 
tion llVe ras &: A rmitaed l2005l . l2006l : [Raymond et al.l 120101 : 


iMatsumura et ah 2013h . Other potential sources of such 


high eccentricities are stellar flybys and Galactic tides; af¬ 
ter a brief period of stellar mass loss, flybys and tides act in 
WD systems in a similar manner to m ain sequence planet ary 
systems, including the Solar system JVeras et al.ll2014rll . 

Because the timescales for orbital perturbations due to 
sublimation are short, and are often strongest at the first 
pericentre passage, they would likely domi nate perturba¬ 
tions from WD radiation. IVeras et al.l (l2015bli demonstrated 
that the ever-dimming WD radiation will shrink and circu¬ 
larize the orbits of sub-metre-sized particles eventually due 
to Poynting-Robertson drag. For the minor planets con¬ 
sidered here, which are km-sized, the dominant radiativ e 
effect would be the Yarkovsky effect dVeras et al.l l2015al 'l. 
When activated, the Yarkovsky effect dominates Poynting- 
Robertson drag, but the timescale for Yarkovsky to act is 
still much longer than that from sublimation and GR. Gon- 
sequently, for volatile-rich comets, sublimation likely alters 
the orbits before Yarkovsky becomes important. In these 
cases, sublimative forces in effect set the initial orbital con¬ 
ditions for other radiatively-induced perturbations. Alterna¬ 
tively, for volatile-depleted minor planets, only the relative 
strength of GR and the Yarkovsky effect would be impor¬ 
tant. 

A sublimative perturbation on the semimajor axis 
of even just a few hundredths of an au may be signif¬ 
icant. In WD systems that contain planets, that differ¬ 
ence can break or initiate a mean motion resonance with 
a planet. Exo-Oort cloud comets which ap p roach WDs 
dAlcock et al.l Il986l: iParriott fc AlcockI Il994 IVeras et al.l 
l2014d:IStone et al.ll2015h may be ejected into the interstellar 
medium. In WD systems which are in the process of forming 
rings and discs from the tidal disruption of other minor plan¬ 
ets, a change in semimajor axis will alter interactions with 
the resulting debris field during each pericentre passage. A 
significant debris field may already exist anyway from the 
end of the giant branch phases of evolution due to YORP 
break-up of km-siz ed asteroids, potent ially extending to tens 
or hundreds of au (IVeras et al.ll2014al . 

Finally, we note that the line between minor planets and 
comets may be fuzzy. Even in the Solar system, the distinc- 
ti on is not always ea sy to mak e (see, e.g., the Introduction 
of ijewitt et al.l [20151 !. In fact, IWeissman fc LevisonI (Il997fl 
claimed that about 1 per cent of the Solar system Oort cloud 


© XXXX RAS, MNRAS OOO.ITlITil 


















































Sublimation of planetesimals in white dwarf systems 


11 


should comprise asteroids. IShannon et al.l ll2015li more re¬ 
cently suggested that this population of asteroids should be 
about eight billion. 


6 CONCLUSIONS 

Hard observational evidence of volatile-rich circumstellar 
material in WD systems motivates studies about how ac¬ 
tive asteroids and comets self-perturb their orbits due to the 
release of volatiles through sublimation and outgassing dur¬ 
ing close pericentre passages with the star. In this paper, we 
have derived the complete unaveraged equations of motion 
dlZliSl and lA2IA6l) in orbital elements due only to sublima¬ 
tion, under some appropriate assumptions. These equations 
are widely applicable to general extrasolar systems assuming 
that planetary and external perturbations are negligible. For 
Solar-system-like sublimation rates, the effect of GR must 
be included (Fig. [!}. We also considered the relative contri¬ 
bution of outgassing in the transverse and nonplanar direc¬ 
tions, and presented corresponding equations (I39I45II that 
also may be applied to exosystems. 

These results, when applied to polluted WD systems, 
ultimately reveal that major planets must exist in order 
to perturb active asteroids into the WD disruption radius. 
Without major planets, these highly-eccentric minor planets 
cannot change their pericentre distances significantly enough 
through the release of volatiles. Inspection alone of equa¬ 
tions (1221, (12H1), (E3) and (I54II reveals why; the semimajor 
axis would change by many orders of magnitude more than 
the orbital pericentre in every case. 
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Table Al. Powers of (1 — e) in the leading order term for the stated expressions expanded about high eccentricity. For W ^ 2, the 
pericentre distance is more resistant to change than any other variable. 



1 

0 

1/2 

1/2 

1/2 

2 

-1 

0 

1/2 

0 

3 

-2 

-1 

0 

-1 

4 

-3 

-2 

-1 
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5 

-4 

-3 

-2 

-3 

6 

-5 

-4 

-3 

-4 

7 

-6 

-5 

-4 

-5 


APPENDIX A: SUBLIMATION WITH DIFFERENT POWER LAWS 


In order to assess how strongly dependent our sublimation findings are on the specific power-law prescription used in equation 
we now rederive the equations of motion under more general considerations. Our physical motivation for the power-law 
exponent was observations from within the Solar system. However, to allow for unknown variations that can arise in the 
populations of small bodies in exosystems, we assume the equation of motion for a sublimating small body is instead given by 


^ sub 


G'(MwD + Mc)r OL /ro\'^r 

r® Me \ r / r 


(Al) 


where W is a, positive integer and a is a constant. This integer requirement will allow us to obtain averaged quantities 
analytically, and further bound the observationally-motivated value with W = 2 and W = 3. This exercise will allow one to 
easily derive the equations of motion for an arbitrary power-law perturbation to the two-body problem algebraically. 

We find that the terms containing the i, Q and the C variables remain unchanged such that 


/ da 
\ dt 


sub 


/ de 
\ dt 


sub 


2arV (l + ecosf)'^ 
Msno'^(l-e 2 )'^+i 
arV (l-|-ecos f)^~^ 


2 Mcna'^+i (1 - e^y 


Cl (cos i cos Q — cos i sin H -|- sin i) — C 2 (cos + sin H) 

Ce (cos i cos H — cos i sin H -|- sin i) — Cs (cos Q + sin H) 


(A2) 

(A3) 


f di\ arV (i-+ e cos f)^ ^ cos (f + uj) . 

— = ———^- - sm i (sm H - cos H + cot i). 

VrfVsub Mcna'^+i(l-e2)'^-5 


(l-e 2 ) 

d^ \ _ arV {1 +e cos f)^ ^sm{f + u}) 


sub Mona'^+i(l-e2)'^-5 

—arV (l + ecosf)^ 


dt J 

dio \ 

/ sub 2Mcnea'^+5 (1 — e^) 

dt 


(sin H — cos Q, + cot i ), 


w 


Cs (cos i cos H — cos i sin fl)+Cg sini+Cr (cos H -I- sin f2)-|-2e sin (/ -|- u) cos i cot i 


2 aro (1 -|- ecos f) 


— (Cip — C2p) , 


de\ 


sub Msno'^(l-e2)'^+5 

P W / 1 I —1 

I = QPq (1+ecos/) 


'^^/sub 2Mcna'^+i (1 — e^)'^ 2 


dzu 

dt 


aro (1 -I- ecos /) 


w-i 


sub 2Mcnea'^'^2(^l — e^)^ 2 


Y {Ctp -I- Cgp) , 


(A4) 

(AS) 

(A6) 

(A7) 

(A8) 

(A9) 


Averaging the above equations presents difficulties analytically and would typically require the use of computer algebra 
so ftware . Even with that tool, we could obtain only the averaged planar equations of motion. We accomplish this task as 
in I Veras! ll2014al ) by computing integrals symbolically where the integrand is a function of an integral power of (1 + ecos/) 
multiplied by single powers of sin (fc/) and cos (fc/), where k is an integer greater than or equal to 1. 

The resulting expressions, denoted by (), are long and we do not present them. Instead, we expand the expressions in a 
Taylor series about e = 1 — e, e -C 1, and consider the leading term in the expansion (denoted by a wide tilde). We hence 
report the lowest-order power of (1 — e) appearing in this expansion in Table ET] 

Because this power is higher for the pericentre distance than for the evolution of any other variable for IT ^ 2, we 
conclude that our main finding is robust for a wide variety of steeper potential sublimation power-law profiles. 
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